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ABSTRACT 



Context. Young planets embedded in their protoplanetary disk interact gravitationally with it leading to energy and angular momen- 
tum exchange. This interaction determines the evolution of the planet through changes to the orbital parameters. 
Aims. We investigate changes in the orbital elements of a 20 Earth-mass planet due to the torques from the disk. We focus on the 
non-linear evolution of initially non-vanishing eccentricity, e, and/or inclination, /. 

Methods. We treat the disk as a two- or three-dimensional viscous fluid and perform hydrodynamical simulations using finite dilfer- 
ence methods. The planetary orbit is updated according to the gravitational torque exerted by the disk. We monitor the time evolution 
of the orbital elements of the planet. 

Results. We find rapid exponential decay of the planet orbital eccentricity and inclination for small initial values of e and i, in agree- 
ment with linear theory. For larger values of e > 0.1 the decay time increases and the decay rate scales as e oc e"^, consistent with 
existing theoretical models. For large inclinations (/ > 6°) the inclination decay rate shows an identical scaling di/dt oc t^. We find an 
interesting dependence of the migration on the eccentricity. In a disk with aspect ratio H/r = 0.05 the migration rate is enhanced for 
small non-zero eccentricities (e < 0.1), while for larger values we see a significant reduction by a factor of ~ 4. We find no indication 
for a reversal of the migration for large e, although the torque experienced by the planet becomes positive when e ^ 0.3. This inward 
migration is caused by the persisting energy loss of the planet. 

Conclusions. For non gap forming planets, eccentricity and inclination damping occurs on a time scale that is very much shorter than 
the migration time scale. The results of non linear hydrodynamic simulations are in very good agreement with linear theory for values 
of e and / for which the theory is applicable (i.e. e and / < H/r). 

Key words, accretion disks - planet formation - hydrodynamics 



1. Introduction 

In the early stages of their formation protoplanets are still 
embedded in the disk from which they formed. Not only 
will the protoplanet accrete material from the disk, it will 
also interact gravitationally with it. Planet-disk interaction is 
an important aspect in the process of planet formation, and 
has been studied already before the dis covery of extrasolar 
planets through l i near analysis ( Goldrei ch & Tremaind Il980t 
iPapaloizou & Lid 119841; IWardI 11986.) . In more recent times 
this was followed up by non- linear numerical simulations us- 
ing a two-di mensional setup ("Brvden et alj 119991: 1 Kiev! 1 19991: 
iLubow et al] |1999; D'Angelo et al. 2002ir a nd later in three 
dimensions as we ll ^'Angelo et alj 120031: iBate et alj 120031: 
ISchaferetalJIIOOl . 

Since detection of the first extrasolar planet in 1995 there 
have been about 200 additional discoveries (for an up-to-date 
list see e.g. http://exoplanet.eu/by J.Schneider). One of the sur- 
prising differences between our own Solar System and these ex- 
trasolar planets are their orbital properties. In the Solar System 
the larger (giant) planets move on almost circular orbits, whereas 
in contrast most extras olar planets move o n eccentric orbits. For 
a recent overview see IMarcv et alj (l2005h . The average orbital 
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eccentricity of observed extrasolar planets is e ^ 0.3. Not only 
are massive planets of several Jupiter masses found in eccentric 
orbits, but so are lighter planets in the mass range where they 
cannot open a clean gap in the disk (M < My„p). 

One possible origin of the eccentricity is through interaction 
with t he protoplanetary d is k. In particular, Goldreich & Sar^ 
(l2003h . ISari & GoldreichI (|2004 estimate that eccentric 
Lindblad resonances can cause eccentricity growth for gap 
forming pl anets. On th e other hand, linear stud i es, fo r ex- 
ample by Ward (1988), Papaloizo u & LarwoodI ( l2000 b and 
Tanaka & Ward (2004), predict eccentricity damping for low 
mass embedded planets. Numerical hydrodynamical simula- 
tions of embedded pl anets tend to show ecc entricity damping 
as well. For example IPapaloizou et al.l (1200 ll) find eccentricity 
growth for bodies on initially circular orbits, but only for 
co mpanions with a mass greater than 10 Mjup. A recent study 
by ID' Angelo et al.1 (l2006l) . however, suggests that modest disk 
induced eccentricity growth can occur for planet masses in 
the Jovian mass range. The origin of the on-average higher 
eccentricities of the extrasolar planets remains to be definitively 
addressed, although recent work indicates good agreement 
between planetary scatte ring simulations and observational data 
(lJuric & Tremainell2007h . 

The linear esti mates of the eccentricity evolut i on of 
embedded plaiiets (lArtymowiczl 119931: IWard & Hahiil 11994 
iTanaka&W^ 12004 " concentrate on small eccentricities and 
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predict exponential decay on short timescales T(,„ ~ (H/r)^T„ng, 
where H/r is the aspect ratio of the disk and r„„g and Tjcc 
the migration and eccentricity damping time scale, respectively. 
iPapaloizou & LarwoodI (|2000|) have also considered larger val- 
ues for e and they find an extended eccentricity damping time 
scale such that de/dt cc e^^ if e > l.lH/r. Additionally they find 
positive torques acting on the planet for large values of e and in- 
terpret this as outward migration. Recently, Cresswell & Nelson 
( I2OO6I) have performed hydrodynamical simulations of embed- 
de d small mass planets and find g ood agreement with the work 
by IPapaloizou & LarwoodI (I2000h . 

The influence of the disk on the incl ination of the planet ha s 
only been analysed in linear studies by iTanaka & Ward! (l2004t) . 
They find exponential damping for any non-vanishing inclina- 
tions on similar timescales as the eccentricity damping. In both 
cases, their results are formally valid only for e, ; <K H/r, but 
numerical experiments jCresswell & Nelsonl l2006') suggest that 
at least their eccentricity damping estimates may be valid for a 
larger range of values and for a variety of sub-gap opening planet 
masses. 

Here we investigate the dynamical influence the protoplan- 
etary disk has on a protoplanet on an eccentric and/or inclined 
orbit. To model the disk dynamics we use a fully non-linear 
hydrodynamical description in both two and three dimensions. 
The orbit of the embedded planet is allowed to evolve due to 
the torques from the disk. While in Kley & Dirksen (2006) the 
possibility of eccentricity growth within the disk due to high 
mass planets (> 1 Jupiter mass) has been studied, we now anal- 
yse the planet-disk interaction of a low-mass planet of 20 Earth 
m asses. Our w ork extends the recent two-dimensional analysis 
of lCressweU & Nelson (2006). We vary the initial orbit param- 
eters and analyse the influence of a nonzero eccentricity and in- 
clination on the migration rate of the planet. In particular we in- 
vestigate if there are important differences between the circular, 
the low eccentricity and the high eccentricity regimes. We anal- 
yse the time scale of the eccentricity and inclination damping 
and compare to previous linear analysis. We also compare the 
migration rates for the different orbital parameters with the cir- 
cular coplanar case and the analytical values from lTanaka et alj 
(I2OO2I) . and determine an experimental upper bound on the valid- 
ity of the exponential decay of ecc entricity and incl i nation mod- 
elled by the analytical estimates of iTanaka & Ward! (|2004|) . 

This paper is organised as follows. In section |2| we discuss 
the disk models and numerical methods. In section|3]we discuss 
the results for planets on initially eccentric orbits located in the 
disk midplane. In sectior|4|we discuss the results for planets on 
inclined and/or eccentric orbits. Finally we discuss our results 
and draw conclusions in section |5] 

2. The Hydrodynamical Model 

To study the evolution of an embedded planet in a protoplane- 
tary disk we treat the disk as a viscous fluid. We consider both 
2- and 3-dimensional (2D, 3D) models which enables us to anal- 
yse the influence of dimensionality directly. Finding the right 
setup for the 2D models to agree with the more complex 3D re- 
sults will enable us to reduce computational efforts in the future. 
The origin of the coordinate system typically is star centred (in 
particular for the 2D simulations and those 3D models presented 
in section nil, but for some simulations it is located at the centre 
of mass of the planet and star system. The z = plane defines 
the disk midplane and the planetary inclination is measured with 
respect to this plane. Many simulations were performed in a ref- 
erence system that would be initially corotating with the planet if 



it were on a circular orbit. This rotation rate is kept constant dur- 
ing the simulations so that we do not have to consider extra ac- 
celerations. Simulations presented in section |4| were performed 
in the inertial frame. 

For the 2D models we use cylindrical coordinates {r,(p,z) 
and consider a vertically averaged, infinitesimally thin disk lo- 
cated at z = 0. The basic hydrodynamic equations (mass and 
momentum conservation) describing the time evolution of such 
a two-dimensional (r, (p) disk with embedded p lanets have been 
stated frequently and are not repeated here (see iKlevll 19991) . The 
2D models presented here are calculated basi c ally in the s ame 
manner as those described previously in Kley ("1998', 1999') us- 
ing the code RH2D. The reader is referred to those papers for 
details on the computational aspects of this type of simulations. 
Other similar models, following explic itly the m otion of single 
planets in disks, hav e been presented bv lNelsori et al. (2000) and 
iBrvden et alj(l2000l) . 

For the 3-dimensional models we work in spherical coordi- 
nates (r, (p, 9) where in the third dimension we use a Gaussian 
density profile, consistent with a locally isothermal equation of 
state. We also set H/r = const. In the vertical direction the grid 
extends over sev eral scale heights. This approach is the same as 
that described in lKlev et all dioOlh and lD'Angelo et all (l2003h . 

2.1. Initial Setup 

The 3D (r, tp, 9) computational domain consists of a complete 
ring of the protoplanetary disk centred on the star. The radial 
extent of the computational domain (ranging from r^m to r^ax) 
is chosen such that there is sufficient space on both sides of the 
planet to reduce possible influence of the boundaries. Typically, 
we assume - 0.4 and r^ax - 2.5 in units where the planet is 
located initially at r = 1 . In the azimuthal direction for a com- 
plete annulus we have ip^in = to (/Jmax = 27r. In the meridional 
direction we take a wedge spanning several disk scale heights: 
Gmin =75° and 0,„a v = 105°. To save computer time, the models 
without inclination were run in half a disk with 9„ax - 90°, using 
symmetry around the midplane. Finally in the 2D simulations we 
use the same radial extent as in the 3D case. 

The initial structure of the disk (density, temperature, and 
velocity) is axisymmetric. For the density p(r, if, 9) we assume 
a Gaussian stratification in the vertical direction with a scale 
length H, and a power law for the radial dependence 

p(r, if, 9) = po(r) exp [-l/2(r9/Hf] (1) 

where we choose poir) cc r"'^ such that the initial surface den- 
sity cr(r) - J pdz falls off radially as a power law with in- 
dex -0.5. If there is no planet, this equilibrium density structure 
will be preserved for constant kine matic viscosity and closed ra- 
dial boundaries (' Kiev et aLll200ll) . The initial velocity is pure 
Keplerian rotation (ur = 0, = (GMt/r)^^^, ug - 0). The tem- 
perature stratification is always given by T(r) oc r"' which fol- 
lows from an assumed constant vertical height H/r - const. For 
these locally isothermal models the temperature profile remains 
fixed and no energy equation is required. We use a constant kine- 
matic viscosity coefficient v. 

To preserve a constant mass in the disk we use reflecting 
boundary conditions at the inner and outer boundary. For test- 
ing purpo ses we apply also dampin g boundary conditions as de- 
scribed in lde Val-Borro et aLl (|2006|) . In the azimuthal direction, 
periodic boundary conditions for all variables are imposed. At 
the upper and lower meridional boundary we also use reflecting 
boundary conditions to preserve the total mass. In those models 
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where only the upper part of the disk is considered we use sym- 
metric boundary conditions at the midplane to simulate a lower 
half of the disk that is symmetric to the upper half. 

2.2. Model parameters 

The computational domain is covered by 131 x 388 x 40 
(A^,. X X Ng) grid cells which are spaced equidistant in ra- 
dius, azimuth and polar angle. Half the number of cells are used 
in the meridional direction when symmetry about the midplane 
is assumed. 

The mass of the planet relative to the mass of the star in the 
different models is 6x10"^, which corresponds to a 20 planet 
for a Solar mass star. We allowed the planetary orbit to evolve 
as a result of the torque exerted by the disk. In these models we 
assume a disk mass inside the computational domain (r = 2.08 
to 13 AU for a Solar type star) of 7.0 X 1Q-^Mq. 

A constant dimensionless kinematic viscosity v - 10"^ is 
used for all models. This corresponds to or = 0.004 for a planet 
at a = 5.2 AU, a typical value for a protoplanetary disk. We set 
H/r = 0.05. 

We let the planet evolve to observe the orbital evolution of 
the system. The motion of the planet is integrated using a fourth 
order Runge-Kutta integrator where the time step size is given 
by the hydrodynamical time step. As a test we have run pure 
2 and 3 -body problems of one star with one or two planets un- 
der the same conditions as in the full hydrodynamical evolution 
(i.e. identical timestep size) and find that the relative total energy 
loss is less than 6 x 10"^ for x 2,500,000 time steps, which is 
equivalent to several thousand periods of the planet. The forces 
from the disk material are calculated in first order and are added 
when the planet positions are updated. Disk self-gravity is not 
included. 

In these models we did not take into account torques (forces) 
from disk material that lies closer to the planet than rtorg - 0.8 - 
1 .0 Rhui where the Hill-radius is given by 



R 



Hill 



(3) 



(2) 



where Op is the actual semi-major axis of the planet orbit. For 
the calculations in section 13.41 we use a smoothed transition for 
the torque cutoff. 

2.3. Numerical issues 

We use two different codes, NIRVANA and RH2D for our cal- 
culations. The numerical method used for both utilises a stag- 
gered mesh, spatially second order finite difference method 
based, w here advection is based on the monotonic transport al- 
gorithm (I van Leeiill977l) . Due to operator-splitting the code is 
semi-second orde r in time . Deta ils of the NIRVANA code have 
been described in IZieded (Il998h . where both the London and 
Tubingen groups have added their own improvements to the 
code. Application of the two dimensional code RH2D to the em- 
bedded planet problem is described in iKlevI (11998). 

The use of a rotating coordinate system in most of our runs 
requires special treatment of the Corio lis terms to ensure angular 
momentum conservation (' KievlfTggsh . 

In calculating the gravitational potential of the planet we use 
a smoothed potential of the form 



- - 



Gm„ 



where s is the distance from the planet. For the smoothing length 
e of the potential we choose the diagonal length of one grid cell 
in three dimensions (this corresponds to 80 % of the Hill sphere 
radius), and to simulate the three-dimensional environment as 
well as possible in two dimensions we use e = Q.6H. 

The viscous terms, including all necessary tensor compo- 
nents, are treated explicitly. 

An important issue in complex numerical hydrodynamical 
simulations is the study of numerical convergence and consis- 
tency. While consistency of our difference equations is satisfied 
through derivation of them from Taylor-expansions, the question 
of convergence is typically addressed through resolution studies. 
As described below, the main results of our simulations are ob- 
tained by running two- and three-dimensional disk-planet sim- 
ulations over several hundreds of orbital periods. Since it is not 
possible to perform resolution studies on all physical setups we 
present in the next sections the results obtained with our stan- 
dard resolution (131 x 388 x 40 for a 3D-setup with inclined 
planet). Additional resolution studies performed on a represen- 
tative sample are described in more detail in the Appendix, and 
briefly at the relevant description of the results in the text. 



3. Eccentricity evolution 

For the study of eccentricity damping of a planet embedded in a 
protoplanetary disk we use here non-inclined orbits and distin- 
guish two different regimes. First we will investigate the regime 
of low non-zero eccentricity and then we will look at high ec- 
centricities similar to those observed in exoplanetary systems. 
The time variation e(t) of the planetary eccentricity is shown for 
different initial eccentricities eo in Fig.[T] 



3.1. Low initial eccentricity 

If we start the planet with a sufficiently low eccentricity (eo < 
0.1) we observe an exponential decay of the orbital eccentricity 
of the planet. 

We find that for the low eccentricity models with eo ^ 0.10 
the decay time t^cc - k/^l is approximat ely 44 - 50 orbits . Using 
linear analysis for small eccentricities iTanaka & WardI (|2004|) 
find that the mean eccentricity change (averaged over one plan- 
etary orbit) is given by 



de/dt _ _ 0.780 

^ ^wave 

with characteristic time 



= 1 



(T„a 



-I 



Cs 



(4) 



(5) 



(,2 + ^2)1/2 



(3) 



where q is the mass ratio between the planet and the star, and o-p 
the local surface density at the planetary orbit. For a planet of 20 
Earth masses at 5.2 AU in a MMSN nebula (o-p - 1490 kg m"^), 
as in our model, the characteristic time is t^^.a^,e = 438 yrs = 37 
orbits. This gives an eccentricity damping time scale of about 
Tecc = twave/O.lE = 47 orbits. This expected exponential eccen- 
tricity evolution is indicated in Fig. [T](upper panel) by the dashed 
curves for eo = 0.1 and (shifted) for eo = 0.2 and 0.3 matched to 
the point where the curves cross e = 0. 1 . In the lower panel of 
Fig. [1] we overlay the eccentricity evolution for different eo such 
that they overlap at e = 0. 1 . The exponential decay phase is very 
similar for the different initial eccentricities (see also below). As 
seen from the plot our calculated e-damping time scale is only 
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Fig. 1. Eccentricity evolution of a 20 niE planet in a 3D pro- 
toplanetary disk as a function of time for different initial ec- 
centricities. We find exponential decay in the low-eccentricity 
regime eo < 0.10, and slower damping for higher eccentrici- 
ties. Top: The dashed lines indicate the exponential linear resul t 
(Tfcc = 47) for this density as derived in lTanaka & WardI (|2004|) . 
Bottom: Time-shifted results for the late exponential phase of 
the simulations with large initial eo ^ 0.1. The dashed curve 
indicates here an exponential decay with Tecc = 44. 



differs slightly, being about 44 orbital periods for e = 0.1 and in- 
creasing later as the eccentricity damps further. This increase in 
damping time relative to that obtained bv lTanaka & WardI (|2004|) 
is probably due to the influence of the gravitational softening in 
our 3D simulations. 

Overall, the linear estimates of iTanaka & WardI (|2004|) are 
therefore in good agreement with our numerical results for ec- 
centricities e < 0.1. Formally it is expected that agreement 
should hold for e < H/r, and for our disks Hjr - 0.05, so the 
linear estimates appear to be good over twice their formal range 
of applicability. 

As can be seen in Fig.|2]an eccentric orbit will also change 
the migration rate of the planet. The nearly straight solid line de- 
notes the migration for a planet on a circular orbit with eo - 0.0. 
Notice that in this case we always maintain the zero eccentric- 
ity. For the eccentric planets the migration rate shows sinusoidal 
variations that are caused by the orbital variations of the torques 
acting on the planet. 
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Fig. 2. Semi-major axis evolution of the planet in a 3D proto- 
planetary disk as a function of time for different initial eccen- 
tricities. For small non-zero eccentricities the migration rate in- 
creases above the e = case. For large eccentricities the migra- 
tion rate is significantly slower. 



An interesting effect can be seen when one looks at the mi- 
gration time scale of the planet. For small initial eccentricities in 
the exponential damping regime, here eo - 0.05 and eo = 0.10, 
we find a slightly faster migration rate for the planet than in the 
circular case. This result is obtained for all codes used in this 
study, and is found in both 2D and 3D runs. Its long term in- 
fluence on the migration rates of protoplanets is not significant, 
however, as the eccentricity damping time scale is shorter than 
the migration time by a factor on the order of (H/r)^. 
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Fig. 3. The models with eo - 0.2 and 0.3 are fitted to a theoretical 
model with e oc e"^. As can be seen this scaling works extremely 
well as long as the eccentricity is larger than about 2.5 H/r. 



3.2. High initial eccentricity 

If we start the planet with a high eccentricity (eo > 0.1) we also 
observe a decay of the orbital eccentricity of the planet, but it is 
slower and the decay rate is no longer exponential, as can be seen 
in the two upper curves in Fig.[T] In fact it fits well with the theo- 
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retical model described in lPapaloizou &Larwoodl ( 120001) . which 
predicts e oc e'^. In Fig.[3]the calculations with eo - 0.2 and 0.3 
are fitted with such a model. The fits are extremely good. The 
characteristic damping time at f = for this non-exponential 
damping with e - Ke'^ is given by 



1.01 



/=() 



1*1 
K 



(6) 



From Fig. |3]we obtain - 1\ for eo = 0.2, and for eo = 0.3 
we find ~ 183 orbital periods, more than a four-fold increase 
over the exponential decay. Also notice that when the eccentric- 
ity falls beneath the value for which exponential damping occurs, 
the lines in Fig. [1] (lower panel) look similar to shifted copies of 
each other, indicating that the planet has no memory of its pre- 
viously higher eccentricity. In the large eo = 0.3 case the addi- 
tional evolution in semi-major axis may account for the observed 
(slight) difference in decay time. 

For these high initial eccentricities we find additionally a re- 
duced migration rate (Fig . |2]i . In Fig . |4]the relative migration rate 
compared to the e = model is plotted. For small eccentricities 
the migration rate may increase by as much as 60% before drop- 
ping off in the manner expected for larger eccentricities. For a 
large eccentricity of e = 0.3 the migration rate is 5 times smaller 
than for a planet on a circular orbit. 
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Fig. 4. The relative initial migration rate as a function of eccen- 
tricity for a planet with mass ratio q = 6x 10^^. For low eccen- 
tricities the migration rate can increase by up to 60%, for high 
eccentricities we find a substantial slowing of the migration. 



In the long term evolution of the migration rate we see in 
Fig. |5] that all models converge to the circular migration rate 
when approaching e = 0, as expected. This turnover to the cir- 
cular migration rate occurs at a residual eccentricity of about 
e = 0.05 ^ H/r, i.e. for the model with eo - 0.10 at around 
t = 20 and for the model with eo - 0.2 at around t - 80. For 
the e - 0.3 model the eccentricity has not dropped enough to 
reach the turnover from the slow migration regime (large e) to 
the standard migration regime (small nonzero e), but will even- 
tually reach the same migration rate as the other models in Fig.[T] 
(see also the two-dimensional models in section l331 l. 

The turnover from the slower migration rate for larger ini- 
tial eo to the circular rate occurs through a brief phase of rapid 
migration when the actual eccentricities are in a range of 0.1 to 
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Fig. 5. The long term evolution of the semi-major axis for dif- 
ferent initial eccentricities for a 20 planet for 3D models. All 
models except the one with the largest initial eccentricity con- 
verge to the same migration rate, and this large e run eventually 
converges also after a longer run time. 



0.05. This occurs for the eo 
to t = 80. 



: 0.20 model during the time f = 60 



3.3. Comparison witli two-dimensional models 

To investigate the influence of dimensionality, we have rerun 
all the previously described models with a low mass eccentric 
planet in a two-dimensional (2D) disk as well. Our initial con- 
ditions for the disk are identical to the models in three dimen- 
sions, only vertically averaged. The grid structure is the one 
described in section 12.21 The evolution of the semi-major axis 
and eccentricity of a 20 Earth-mass planet starting on an ec- 
centric orbit with initial eccentricity between and 0.3 for the 
two-dimensional models is shown in Fig. |6] In the top panel 
the semi-major axis is shown. When the potential smoothing de- 
scribed in section [23] takes the value e = 0.6//, the migration 
rate in the two-dimensional disk agrees well with for the migra- 
tion rate in three dimensions obtained using linear theory. The 
da shed referenc e lines refer to these linear 3D values as given 
bv lTanaka et al.l (l2002h . The obtained migration rate agrees with 
our 3D-results as shown above. Similarly to the 3D case, the 
models with small eccentricities e < 0.10 migrate initially at a 
faster rate, while those for larger e are migrating at a slower rate. 
Another feature that is preserved from the three-dimensional cal- 
culations is the sharp change in migration rate when the eccen- 
tricity drops below x 0.10. For the eo = 0.30 model this happens 
after approximately 240 orbits, for the eo = 0.20 model after ap- 
proximately 70 orbits. As described above, at this point during 
the evolution the eccentricity damping changes from the slower 
e oc e"^ behaviour to the faster exponential behaviour during 
which migration is slightly accelerated. Only when the eccen- 
tricities become very small the migration is slowed again and 
the standard (circular) rate is approached for all models. 

There are also some distinct differences. In the bottom panel 
of Fig. |6] the eccentricity is plotted as a function of time. 
Over-plotted (dashed curves) are the linear damping rates by 
Tanaka & Ward (2004), where our damping is again slightly 
faster. Our damping rates for the 2D and 3D cases agree very 
well for the smaller eccentricities while for the large value eo = 
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As a next step we analysed the role of potential softening 
in more detail. While the use of a larger softening (e - 0.6//) 
to approximate 3D disk structure is a well-known approach in 
2D modelling, for sufficiently eccentric orbits in which a large 
region of the disk is sampled, our results indicate that a simple 
prescription for potential softening may no longer be sufficient 
to obtain agreement between 2D and 3D simulations. This may 
be due in part to the fact that for highly eccentric orbits the planet 
is crossing the resonances that drive the disk-planet interaction, 
such th at a strong sensitivity to the so ftening prescription is ex- 
pected (iPapaloizou & Larwoodll2000l) . Further calculations that 
we have performed indicate that a considerably more complex 
prescription for softening the gravitational potential will be re- 
quired if 2D models are to be constructed that agree with 3D 
models for both high and low orbital eccentricities. Clearly, a 
softening that is spatially fixed and does not vary with radius 
(which H does in fact) will not be correct, a conclusion which 
applies to the 3D case as well in case the simulation is resolu- 
tion limited. In our case at r = 1 the standard grid resolution 
yields a diagonal length of A « QSRhui for one gridcell. The 
torque cutoff r,orq should also scale with the radial distance of 
the planet from the star and not be a function of semi-major axis 
alone. We would like to point out that increasing the resolution 
of our models in the two- and three-dimensional case changes 
the time evolution of the semi-major axis only marginal, but lead 
to a slightly faster eccentricity damping. However, the observed 
discrepancy in timescales between the 2D and 3D case is not al- 
tered. See Appendix for a more detailed discussion on resolution 
issues. 



3.4. Dynamics of the How 



Fig. 6. Top: Semi-major axis as a function of time for q - 
6 X 10"^ in a 2D disk for different initial eccentricities. For 
very high eccentricities the speed of migration is significantly 
reduced as in the 3D case. Dashed lines indicate the result from 
linear analysis for 3D disks. Bottom: Eccentricity as a function 
of time for the same runs. The das hed lines indicate the linear 
results from lTanaka & WardI (|2004 with t^cc = 47. Our 3D re- 
sults are over-plotted using the small black symbols. For larger 
initial eccentricities 2D and 3D runs begin to differ; by about 
20% for the e^) - 0.3 case. Possible reasons for this behaviour 
are discussed in the text. 



0.3 the 3D run gives a significantly faster damping. However, 
the generic features and scaling of the damping process for large 
Co is captured correctly in both 2D and 3D runs, and the differ- 
ences in the eccentricity damping are at the 20% level even for 
the largest e^). 

To examine the origin of this discrepancy, we performed ad- 
ditional 2D and 3D models with larger disks, 0.25 < r < 3.5, 
while preserving grid resolution, to rule out the possibility that 
boundary effects are the cause. A small difference in eccentric- 
ity reduction was observed at the transition to linear damping, 
dejdt oc e, with longer non-linear and shorter linear e-folding 
times than for our standard disk model, but the difference is only 
a few percent and not of the order observed in Fig. |6l Runs 
with enhanced damping at the radial boundaries to reduce re- 
flections did not change the results considerably, some damping 
is clearly required however, to avoid unphysical influence due to 
the boundaries. 



To analyse the dynamical structure of the flow we consider the 
change of the flow field in the case of an eccentric planet for a 
2D case. In Fig. |2]we display the density structure of the disk 
with an embedded q - 6x 10"^ planet on an eccentric orbit with 
fixed e = 0.1 at four different phases in the orbit, separated by 
1/4 orbits, ranging from f = 10 to f = 10.75. At f = 10 the planet 
is at apoastron (top left) and at f = 10.5 at periastron (bottom 
left). In contrast to the evolution with zero eccentricity where 
two stationary trailing arms exist, one in the outer disk (r > r^) 
and one in the inner disk (r < r^), in this eccentric case spiral 
arms and additional flow features appear and disappear period- 
ically in phase with the orbit. During periastron a pronounced 
outer spiral attached to the planet is visible while at apoastron 
this has changed to an inner spiral. This periodic shift of the spi- 
ral arm strengths will result in a corresponding periodic variation 
of the torque acting on the planet. We note also that a significant 
density enhancement appears in the close vicinity of the planet. 
At apoastron it clearly lies in front of the planet (and thus ex- 
erts a strong positive torque), and at periastron it lags behind 
the planet (exerting a strong negative torque). This flow feature 
appears to arise because when the orbit is eccentric the flow in 
the planet vicinity becomes similar to a Bondi-Hoyle flow. At 
apoastron the planet moves more slowly through the gas, and so 
is over-taken by the disk matter on trajectories that lie both in- 
side and outside the planet orbit. These flow lines are distorted 
by the gravitational field of the planet and come to a focus in 
front of the planet, forming the high density feature seen in the 
top left panel of Fig. Q The same effect happens in reverse at 
periastron, leading to a high density feature that forms behind 
the planet. These features become dominant in determining the 
torques experienced by the planet. 
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Fig. 7. Density contour plot for a 20 Earth-mass planet on an fixed e = 0. 1 orbit imbedded in a 2D disk at four difl'erent times during 
one orbit. The snapshots are separated by 1/4 orbits. 



This is exemplified in Fig. [8] where we display the variation 
of the torque (see Eq. |9]) acting on the planet and its radial dis- 
tance form the star as a function of time for two different plan- 
etary eccentricities (top: e - 0.1, and bottom: e - 0.3). The 
top plot (e - 0.1) refers directly to Fig. |2l Clearly in apoastron 
(strong inner spiral and leading high density feature) the total 
total torque is positive, while at periastron (strong outer spiral 
and lagging high density feature) the contribution is negative. 
Additionally there appears to be a small phase shift between the 
distance and the torques. Typically, for an embedded protoplanet 
the direction of migration and its magnitude is attributed to the 
sign and value of the torque acting on it. For the highly eccen- 
tric case (e = 0.3) the average torque is clearly positive and one 
might expect outward migration. 



However, this conclusion must be corrected for possible 
changes in the eccentricity of the planet. For a planet on an ec- 
centric orbit its angular momentum Lp is given by 

Lp = mp VCM.fl Vl -e2 (7) 

and the rate of change of the semi-major axis and eccentricity 
can be obtained from 

- i ^ _ i _ ^disk 

Lp 2 fl \ - e Lp 

Here Tdisk is the total torque exerted by the disk onto the planet 
rd,sk = r (rp X F)\ df (9) 

JDisk 
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Fig. 8. The change of the distance of the planet from the star and 
total disk torque acting on the planet as a function of time, where 
the torque is rescaled appropriately. The planet is held fixed on 
its orbit. Top: e = 0.1 and Bottom: e - 0.3. 

where denotes the radius vector from the star to the planet, F 
the (gravitational) force per unit area between the planet and a 
disk element (at location r from the star), and df the surface ele- 
ment. Eq. (ISj implies that a positive torque may also result in ec- 
centricity damping rather than outward migration. From our sim- 
ulations we are able to disentangle the different contributions. 

The evolution of the semi-major axis and eccentricity 
changes are displayed in Fig. |9]for two different initial eccen- 
tricities in the time interval 10 to 13 orbits. After this time the 
flow has equilibrated and the contribution of the individual terms 
in Eq. ([8]l can be analysed. In these calculations we have used a 
smoothed torque cutoff to reduce the noise in the curves. In both 
cases (eo -0.1 and eo - 0.3) the sum of the two terms for a 
and e in Eq. ([8]l (thick dashed line) equals exactly the total disk 
torque as given by the squared symbols. The periodic behaviour 
of all quantities is again clearly visible, they all fluctuate around 
their zero value. Quite clearly, the contribution of the e term has 
comparable magnitude to the a value. Hence, a positive average 
total disk torque as for example in the e - 0.3 case does not 
necessarily imply an outward migration of the planet, since the 
total torque is "shared" between semi-major axis and eccentric- 
ity change (see Eq. [8]). A positive torque implies only that the 
angular momentum of the planet has to increase. For an eccen- 
tric orbit this can be achieved in two ways, either by increasing 
a (outward migration) or by reducing e (circularization). Indeed, 



Fig. 9. The time changes of the semi-major axis (l/2a/a, solid 
curve), the eccentricity (-e^/(l - e^)e/e, thin dashed line), and 
the normalised torque (Tdisk/ip) for an evolving eccentric planet 
in a 2D disk. The sum of the first two contributions is given by 
the symbols. Top: eo = 0.1 and Bottom: eo = 0.3. 



for the semi-major axis to increase the planet's total energy must 
increase. As described below, we find that it in fact decreases 
even though the angular momentum increases through the posi- 
tive torque. 

3.5. Migration rate 

As noticed in the previous sections the migration rate depends on 
the eccentricity of the embedded planet and can vary by about 
60% (cf. Fig. m, an effect seen in the 2D as well as in the 3D 
simulations. For small eccentricities, that are less than or equal 
to about twice the disk aspect ratio (e < 2H/ r), we find an in- 
crease in the migration rate. For larger values of e the rate is 
reduced with respect to the circular case, but always directed in- 
ward. To analyse this effect we have performed additional simu- 
lations in 2D with varying eccentricities. To measure the values 
of the torque, the planet was held again on a fixed orbit. 

In Fig.[TO]we display the the evolution of the torque (top) and 
mechanical power (bottom) as a function of time over 2 orbital 
periods for different eccentricities. Here the energy change per 
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time (power) of the planet due to the work done by the gravita- 
tional forces of the disk is given by 



JDisk 



(10) 



The energy of the planet depends only on the semi-major axis a 
and is given by 



1 GM.mp 



(11) 



2 a 

For the energy loss and semi-major axis change we then obtain 



isk 



Restriction to circular orbits with e = yields 



|£pl " ip 

or for the normalized torque and power 



^disk 2 -^disk 



(12) 



(13) 



(14) 



This relationship which can be directly verified from Fig. [TO] 
where the normalised torque and energy loss is displayed as a 
function of time for various eccentricities. For low eccentricities 
(e = 0: inverted triangles, e - 0.02: open diamonds) the shape 
of the curves are in good agreement, and only scaled by a fac- 
tor of 2. However, for eccentric non-circular motion this simple 
relation is not satisfied anymore. 

Already at very low eccentricities the variations of the power 
exceed the mean value for e = by a large margin. For zero 
eccentricity the value for the power is about -0.09 in the dis- 
played units (solid curve with inverted triangles). For e - 0.02 
(light dashed curve with open diamonds) we find an amplitude of 
0.37, and for e - 0.05 (dark dashed curve with filled diamonds) 
it has increased to 1.00 (bottom plot in Fig.fTOl). Hence, a small 
asymmetry in the power may lead to a substantial change in the 
migration rate. From Fig. [TO]it is clear that in particular during 
periapses {t = 10.5 and 1 1.5) the planets experiences a large en- 
ergy loss which reaches a maximum for an eccentricity around 
e = 0.10. As a consequence the mean value (of Pdisk) has signif- 
icantly dropped below the circular case, leading to the enhanced 
inward migration found previously. For larger e the amplitude 
of the power variation remains approximately at the same value 
but becomes more symmetric, leading to a reduction in the mi- 
gration rate. As the orbit-integrated power remains negative for 
the high eccentricity cases (i.e. e ^ 0.3), the migration remains 
inward even though the orbit integrated torque is positive. 

The variation of the torque for non-vanishing eccentricities 
also increases substantially above the zero eccentricity value 
(Fig. [To] top panel). Here, for larger values of e above about 
e = 0. 1 the mean value becomes clearly positive. However, due 
to the corresponding change in eccentricity and energy the mi- 
gration is still directed inwards (cf. Figs.[8]and|9]). 

To illuminate this effect from a different perspective we have 
performed a set of simulations where the planet was not allowed 
to move and remained on a fixed orbit during the simulations. 
We used 22 different values for the eccentricity ranging from 
e = to e = 0.40 and ran all models up to a final time of 100 
orbits where the torques and the power acting on the planet have 
reached an equilibrium on average. For the last 20 orbits (from 




Fig. 10. Top: Normalised Torque (Tdisk/ip) vs. time for various 
eccentricities, where the orbit is not allowed to evolve. Bottom: 
the corresponding energy loss (Pdisk/|£^pl)- The values have been 
amplified vertically by the same amount as in Fig.|8] 



f = 80 to 100) we calculate the time average of Tdisk and f disk- 
All these models are performed using the two-dimensional setup 
and utilize a 256 x 700 grid with a logarithmic spacing in the 
radial direction. Additionally, to increase performance th e runs 
use the FARGO-algorithm for differentially rotating flows (Massed 
2000)- The results are displayed in Fig.[TT]where the inverted tri- 
angles refer to the torque and the solid dots to the power acting 
on the planet, both given in dimensionless units. Clearly, as al- 
ready found through the previous analysis an increase in the ec- 
centricity leads to a non-monotonic behaviour of both Tdisk and 
^'disk- For circular motion both are identical and negative, leading 
the the well known inward migration of the planet. Upon increas- 
ing the eccentricity both quantities initially drop even more and 
increase later again. The power remains negative for all eccen- 
tricities leading to the inward migration which only slows down 
for larger e. This leads to the behaviour of migration time scales 
as shown in Fig. |4] above, which are shortest for e » 0.08. The 
torque becomes positive for all e > 0.08 and reaches a maximum 
foie^ 0.17. 



3.6. Changing the disk's density profile 

All the previous simulations have been peformed using a sur- 
face density profile which varies as S(r) oc r '', where we 
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Fig. 11. Dimensionless torque (Tdisk) and energy loss (power, 
^'disk) experienced by the planet for different eccentricities. The 
values are obtained by keeping the planet on a fixed orbit and 
using time averages over 20 orbits (from f = 80 to 100). 



used p - 1 12 which is the equilibrium profile for a constant 
kinematic viscosity and closed radial boundaries (M - 0). 
The corotation torques acting on the planet and the resulting 
changes in migration speed depend (for circular orbits) on 
the density slope of the disk (Tanaka et al. 2002; Masset et alJ 
[2006). To test the possible influence of the density slope on 
the orbital evolution of eccentric (non-inclined) planets, we 
have varied the value of p from to 1.5 for our highest value 
of the eccentricity, e - 0.3. The results of our simulations are 
displayed in Fig.[T2j For an increasing density slope the mi- 
gration occurs monotonically at a faster rate, in agreement 
with standard Unear estimates for the torques dTanaka et alj 
[2002*). The eccentricity evolution is much less affected, the 
largest change occurs only for the steepest slope p - 1.5. 
These results indicate that in the case of eccentric plan- 
ets the importance of the corotation torques are greatly re- 
duced. This is probably because during each orbit, the planet 
moves a distance radially that is larger than the width of the 
horseshoe region, where the corotation torque is generated. 
Clearly it is expected that the corotation torque will weaken 
significantly under these conditions. In all four cases with 
different p the average torque is clearly positive while the 
power is negative. As in the previous p - 1/2 case this leads 
to inward migration and eccentricity damping. In particular 
our results concerning the p ositive torque are in agr e ement 
with the earlier findings of 'Pap aloizou & LarwoodI ( l2000l) 
who used p - 1.5, but fixed the planetary orbit and did not 
measure the power. 

4. Inclined orbits 

In addition to the eccentric planetary orbits we now study the 
additional degree of freedom provided by inclined orbits. We 
begin by considering inclined, circular orbits, before discussing 
the general case when both e and i are non- vanishing initially. 

4.1. Circular Orbit 

If we start a low mass planet on a circular but inclined orbit 
the orbital inclination will be damped (cf. Fig. [T3] ). For / ^ e. 
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Fig. 12. Top: Time-change in semi-major axis for different sur- 
face density profils of the disk, given by 2(r) oc r^P. Bottom: 
Corresponding eccentricity evolution. 

the time scale for inclination damping is somewhat longer than 
for the eccentricity damping, but of the same order of magni- 
tude. Again we find two different regimes. For ; < H/r (/ in 
radians) the planet remains inside the main body of the disk 
and the damping is exponential: di/dt oc /. The damping rate 
(averaged over one planetary orbit) for a small planet mass 
and small inclinations is obtained through linear calculations by 
lTanaka&Wc^(l2004l) as 



di/dt _ 1 _ 0.544 

with tn-ave defined in equation (|5]l. For our simulations this gives 
an inclination damping time scale t,„c of about 68 orbits. The 
two lower curves in the upper panel of Fig. [13] show the incli- 
nation damping for small initial inclinations (2° and 5° respec- 
tively). The fitted dashed line for /q = 2° coiTesponds to an ex- 
ponential damping with a timescale of 90 orbits. This obtained 
T,„c is larger than the linear estimate (Tanaka & Ward 2004) by 
the same factor as the eccentricity damping time r^cc- 

For higher inclinations we find a slower inclination damping 
rate as can be seen in Fig. [13] for /(> = 10°. For our disk pa- 
rameters, the damping rate departs from being exponential for 
values of / in the interval 6° < i < 8°. Just as for the damping 
of eccentricities , we th erefore find that the linear calculations of 
iTanaka & WardI (|2004|) provide a good approximation for over 
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Fig. 13. Orbital inclination (top panel) and semi-major axis 
(lower panel) as a function of time for different initial values 
of i for a 20 Earth-mass planet on a circular orbit. We find ex- 
ponential decay of / as long as / < 6°; the thin dashed line su- 
perimposed on the lowest iq - 2° curve corresponds to an expo- 
nential with a damping timescale t,„c - 90. For higher inclina- 
tions the damping timescale increases dramatically. A very good 



fit is found for a model with di/dt 
io = 10°). 



,-2 



(thin dashed line for 



twice their formal range of appUcability, i.e. up to / ~ 2H/r 
- where here / is measured in radians. For / > 8° the damping 
rate strongly deviates from being exponential, and the best fit to 
the model with the highest initial inclination of 10° is given by 
di/dt cc i^^ (upper dashed line in Fig.[T3]). Interestingly, this de- 
pendency is identical to the behaviour of the eccentricity damp- 
ing for large eccentricities. However, here the reduced damping 
comes from the fact that for inclinations significantly larger than 
H/r the planet's contact with the disk is reduced, and its velocity 
relative to local disk gas as it crosses the midplane is increased. 

As with increasingly eccentric planar orbits, the migration 
rate is observed to decrease with increased inclination, as shown 
in Fig. [13] (lower panel). The effect of an increased inclination 
is weaker than for a comparable value of eccentricity: circular, 
inclined orbits lead only to weaker torques on the planet, and 
not the possibility of torque sign reversal. Migration rates return 
smoothly to values close to the planar, circular rate once the in- 
clination drops beneath about 5°. 
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Fig. 14. Semi-major axis, eccentricity and inclination as a func- 
tion of time for a 20 niE planet with different initial eccentricities 
and inclinations. In the top panel the time coordinate has been 
stretched to show effects more clearly. 



4.2. Non-circular Orbits 

Finally we investigate several models in which both the initial 
eccentricity and the initial inclination are nonzero with < eo < 
0.30 and < /o < 8°. Each system was initialised with the planet 
at apocentre above the midplane, with the longitude of pericentre 
o) - jt/2. Our results are summarised in Fig. [14] 
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4.2.1 . Small e and small / 

We begin by considering the evolution of orbits for which the 
planets are initiated to have small values of eccentricity and in- 
chnation (i.e. eo ^ HjR and /q < HjR, where /q is measured in 
radians). 

Examining the top panel of Fig. [14] we can see the migration 
behaviour for runs with (eo = 0.01 and /q = 0.3°) and (eo = 0.05 
and ;o - 4°). It is clear that for small values of both eo and iq 
the migration behaves as if the orbit was essentially uninclined. 
Small values of the inclination such that / < HjR, with / mea- 
sured in radians, cause migration to be slowed, but only very 
slightly. 

Examining the middle panel of Fig. [141 we can see that small 
values of eo and /o lead to exponential decay of the planet eccen- 
tricity, with the rate of damping only slightly increased by the 
small inclination. Examining Fig.[T]we see that the time taken 
for the eo =0.1 case to halve its eccentricity is ^ 35 orbits, and 
the time required for the eo - 0.05, i'o = 4° case shown in Fig. [14] 
to also halve its eccentricity is approximately 40 orbits. 

The lowest panel in Fig. [T4l shows the inclination evolution. 
In common with the eccentricity damping, we find that the in- 
chnation damping for an orbit with small initial values of e and 
; occurs exponentially on a time scale very similar to the case 
when the orbit is circular but inclined. 



4.2.2. Small / and large e 

We now consider the orbital evolution when the inclination is 
small a < H/R) and the eccentricity is large (e > HjR). This 
is represented by the runs shown in Fig. [14] with (eo = 0.3 and 
i'o - 0.5°) and (eo = 0.3 and i'o - 4°). Also shown for reference 
is a case with (eo = 0.3 and /o - 0). The migration is shown in 
the top panel. We have already shown in previous sections that 
an uninclined, high eccentricity orbit undergoes slower inward 
migration than a circular orbit does. We see that increasing the 
inclination increases the migration time further A very small 
inclination (/ = 0.5°) has very little effect, but an inclination of 
; = 4° begins to make a noticeable difference, decreasing the 
migration rate by ^ 25%. 

The eccentricity evolution of these runs is shown in the mid- 
dle panel of Fig. [14] The discussion presented in previous sec- 
tions showed that an uninclined, highly eccentric orbit shows an 
eccentricity decay rate dejcit cc e'^. We can see that a very small 
inclination (/ - 0.5°) has essentially no effect, but an inclina- 
tion of / - 4° extends the damping time by about 25%, while 
preserving the shape of the curve of e versus time. 

The inclination evolution is shown in the lowest panel. If we 
consider the evolution of the runs with /o - 4°, then we see a 
fairly dramatic change in the inclination damping once eo >> 
HjR. For small values of eo we see the usual exponential decline 
of /. When eo is large, however, the inclination damping rate 
is small (even for small inclination) and almost constant with 
superimposed oscillations until the eccentricity falls below e ^ 
0.1, after which it recovers the exponential decay. For example, 
the jQ = 4° and eo = 0.3 case shows an inclination decay time 
scale iolidildt) ~ 700 orbits initially. Once the eccentricity has 
damped to e < 0.1 after approximately 200 orbits, the inclination 
declines exponentially on an e-folding timescale of ^ 90 orbits. 
The period of the oscillations seen in the inclination corresponds 
approximately to the precession period of the nodal line of the 
planet. 



4.2.3. Large / and small e 

We now consider the evolution when the inclination / is large 
and the eccentricity is small, which is represented in Fig. [T4lbv 
the run with (/o = 8° and eo - 0.05). The migration is shown 
in the top panel, and we can see that compared to the run with 
(i'o = 4° and eo = 0.05) the migration rate for the more inclined 
planet is slower by about 40%, such that a doubling of the initial 
inclination leads to a migration rate that is almost halved. After 
just over 220 orbits the migration rate approaches that for un- 
inclined orbits as the inclination has decreased sufficiently after 
this time. 

The eccentricity evolution for this run is shown in the middle 
panel of Fig. [14] Comparing this run with the one with (/o = 4° 
and eo - 0.05) we see that having a large inclination substan- 
tially increases the damping time. Close to the beginning of the 
simulation the low inclination run has its eccentricity reduced 
by about 50% within ^ 40 orbits, whereas the more inclined 
run shows only a 20% change in eccentricity after this time. 
The damping of eccentricity also shows a strong deviation from 
exponential decay for the high inclination run, with the eccen- 
tricity decay rate being almost constant until the inclination has 
damped down to ^ 4°. Once the inclination has declined down 
to this value the eccentricity decay becomes exponential. 

We have already discussed the inclination evolution for cir- 
cular orbits, and have shown that for large values of / the in- 
clination damping rate is no longer exponential, but scales as 
di/dt oc r^. The addition of a small eccentricity (eo = 0.05) 
barely changes the decay rate for large inclinations, as shown in 
the lowest panel of Fig. [14] The effect of this small eccentricity 
is simply to lengthen the damping time by about 10%. 

4.2.4. Large / and large e 

We now consider the evolution of orbits with both large incli- 
nation and eccentricity. The run with (eo - 0.3 and / - 8°), 
plotted in Fig.[T4lillustrates this case. The migration rate, shown 
in the top panel, is seen to be reduced substantially compared to 
both the cases with (eo - 0.3 and /q = 4°), and (eo - 0.05 and 
io - 8°). An increase of inclination from 4° to 8°, with eccen- 
tricity e - 0.3 leads to a ^ 50% reduction in the initial migra- 
tion rate. However, we note that the rapid damping of inclina- 
tion and eccentricity means that this reduction in migration rate 
is short lived, and can at most only increase the total migration 
time through the disk by < 10% i n the absence of a mechani sm 
to maintain e and / at large values dCresswell & Nelsonll2006l) . 

The eccentricity evolution shown in the middle panel demon- 
strates that when eo = 0.3, increasing the inclination from 4° 
to 8° causes a significant lengthening of the initial eccentricity 
damping time (by about a factor of 2). As for the other high ec- 
centricity/inclination runs, the damping is no longer exponential 
while e > 0.1 or / > 4°, but becomes so once the eccentricity 
and inclination fall beneath these values. 

The inclination damping shows similar behaviour to the ec- 
centricity damping for this large e and / run, although the de- 
pendence of di/dt on / is slightly weaker than that seen for 
de/dt. Increasing the inclination from ; = 4° to / = 8° with 
eo - 0.3 decreases the initial inclination damping rate by about 
40%, with little change to the overall shape of the curve. We see 
that both the eccentricity and inclination decay on similar time 
scales (within about 300 orbits for e and 350 orbits for /), and 
once they have reached small values they under go exponential 
decay toward zero. 
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5. Conclusions 

We have performed fully non-linear 2- and 3-dimensional hy- 
drodynamical simulations of a planet embedded in a protoplan- 
etary disk. The planet was allowed to change its orbit due to 
the torques from the disk material. We investigated the orbital 
evolution of the planet in the linear low mass regime (Type I 
migration) for which analytical studies of the eccentricity and 
inclination damping rates are known. 

Concerning the damping of eccentricities and inclinations, 
for small values e < 0.1 and / < 5° we find very good agree- 
ment between our non-linear results and existing linear theory 
dTanaka & Ward 2004). The damping occurs for both, / and e, 
exponentially with a time scale t oc (Hlr)^T„,ig. The quantita- 
tive agreement of our numerical results with the estimates of the 
linear calculations for the eccentricity and inclination damping 
are at the 20-30% level, which is very encouraging consider- 
ing that one always has to use both gravitational smoothing and 
torqu e-cutoff within the planet Hill sphere in numerical simula- 
tions. |P^aloiiou^L^wool (|200^ have shown through linear 
calculations that the absolute magnitude of the torque depends 
on the value of the gravitational smoothing e, and the slightly 
longer damping times that we find in 3D simulations is prob- 
ably an indication that softening is playing a role. In addition, 
iMasset et al.l (12006 ) have suggested that departures from the lin- 
ear results may occur for 20 Earth mass planets that we have 
considered, so we should probably not expect perfect agreement 
with the results of Tanaka & Ward ( 2004.) . 

For our adopted HIr - 0.05, the critical values for e and / (in 
radians) lie approximately at IHjr. Above those critical values 
the time behaviour changes for both, eccentricity and inclina- 
tion. In both cases we find a damping following dyldt oc y^^ 
with y G {e , i], which (for the eccen tricity) has been suggested 
recently by iPapaloizou & Larwoo IdlOOQt). For the inclination 
we attribute this behaviour to the fact that on highly inclined or- 
bits the planet loses contact with the disk and experiences only 
reduced damping, while for the eccentricity it is has been at- 
tributed to th e varying planetary velocit y with respect to the disk 
material (Pap aloizou & Larwoodll2000h . 

For the migration rates as a function of eccentricity we find 
some surprising results: For low eccentricities in the exponential 
damping regime we find faster migration than in the case of zero 
eccentricity - an increase of up to 60% is observed. Through a 
detailed analysis of torque and energy loss balance during the 
planetary orbit we attribute this increased migration rate to an 
increase in the energy loss during periapse of the planet which 
reaches a maximum at e a: 0.1. For larger eccentricities the mi- 
gration rate is reduced significantly below the circular case but 
is still directed inward even for e - 0.3, even though the average 
torques are clearly positive. Migration can still be inward in this 
case because the torque largely goes into damping of the eccen- 
tricity while the planet continues to lose energy on average. If a 
situation were to arise where this energy loss could be compen- 
sated for, for example by gravitational interaction with surround- 
ing planets, then in principle the posi tive disc torque could drive 
outward migration. Simulations by Cresswell & NelsonI (l2006l) 
have explored this possibility, but find that eventually the eccen- 
tricity of all planets damps and inward migration of the planetary 
swarm occurs. 

We compare our 3D results in detail with corresponding 2D 
simulations in the case of vanishing inclinations. For coplanar 
orbits and small eccentricities the results of the 2D and 3D sim- 
ulations are in excellent agreement with respect to the migra- 
tion rate and the eccentricity damping time scale. The increased 













stc 

MedRes 
HighRes 


: 128x38 
: 256x76 
: 384x115 


\ — . — 

J 































































































































50 100 150 200 250 300 350 400 

Time [Orbits] 













std 
MedRes 
HighRes 


: 128x38 
: 256x76 
: 384x1 15 


1 — ■ — 
3 9 
I . 



































































































50 100 150 200 250 300 350 400 

Time [Orbits] 

Fig. 15. Semi-major axis and eccentricity as a function of time 
for a two-dimensional disk with a 20 planet and an initial 
eccentricity eo = 0.3 for different grid resolutions. 



migration for small non-zero eccentricities is also reproduced. 
For eccentricities e > 0.25, agreement between the 2D and 3D 
is poorer. While the qualitative behaviour is still very similar, 
we find differences in the time scales. For example, the eccen- 
tricity damping rate between 2D and 3D differs by 10-20% for 
the eo = 0.3 case. As the planet samples a large range of cell 
sizes and the interaction becomes more sensitive to the smooth- 
ing length selected. This issue can be solved if detailed knowl- 
edge of the problem under examination is known beforehand, 
but with the increased availability of high-performance and par- 
allel computing facilities, it is questionable as to whether the 
results are worth the investment. Hence, the evolution of em- 
bedded planets can be studied in 2D, for eccentricities of up to 
several scaleheights, if the right smoothing length is chosen for 
the potential. For larger values of eccentricity, e-folding times 
are typically overestimated. 

For the combined general case, non-zero eccentricity and in- 
clined orbits, we find that provided the eccentricity or inclination 
remains small, the previous formulae remain a good approxima- 
tion of a planet's behaviour Inclination produces a weaker de- 
parture from the circular, planar case than the same value of ec- 
centricity, and non-linear effects set in sooner for lower e when 
both parameters are non-zero. The migration rate may be mod- 
erately (factors of ~ 2 - 3) reduced if both e and / are towards the 
upper limits of this range. Once both parameters enter the non- 
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Unear regime however, cross-terms become significant in both 
the damping and migration, and the previous correspondences 
(de/dt oc e"2, di/dt oc r^) are lost. Secular effects from the disk 
can lead to moderate, short-term oscillations in e and /. 

Our results indicate that the interaction of a small mass 
planet (here 20 iiie) with its protoplanetary disk always leads to 
a rapid damping of its eccentricity and inclination. Additionally, 
the main effects and timescales are captured very well by a 
two-dimensional approach which greatly simplifies the compu- 
tational effort. Only if additional perturbers, such as other plan- 
ets or a stellar binary companion, are present may the eccentric- 
ity/incUnation of the planet and the disk be excited. 
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Appendix A: Numerical convergence 

To test the issue of numerical convergence we have run a subset 
of the previously described models with a higher grid resolution. 
In all simulations we use a torque cutoff of r,,,^, = 0.8/?//,// and 
(in the 3D cases) a smoothing length of identical extension. 

The first set of simulations refer to the two-dimensonal setup 
(2D-disks). From our standard resolution in 2D with Nr x - 
128 X 384 we first doubled the number of gridpoints in each 
dimenension to a medium resolution of 256 x 768 and then in- 
creased it by again a factor of V2 = 1.41 to 384 x 1152, our 
high resolution case. The results for an initial eccentricity of 
eo = 0.3, as displayed in Fig. [15] indicate first a qualitatively 
identical behaviour for all cases, with only a minor shortening of 
the eccentricity damping timescale with higher resolution, which 
amounts to about 15% difference in between the standard and 
the HiRes case. The differences in the migration can be attributed 
to the changing eccentricity damping, as in the initial phase of 
the evolution the migration rate is nearly identical for all three 
models. By chosing the highest value for the eccentricity we in- 
tended to pick out the most extreme case, and we expect the 
variations to be smaller for lower eccentricities. 

Running a set of three-dimensional simulations on a typ- 
ical parameter set with eo = 0.2, /q = 5° with ob- 
tain the results displayed in Fig. [16] Here, low resolu- 
tion refers to {Nr,Ng,N^) = (132,40,390), medium resolu- 
tion to (Nr,Ng,N^) - (184,56,560) and high resolution to 
(Nr,Ng,N^) = (264, 80, 800), i.e. the three cases differ by about 

a factor of V2 in the number of grid cells in each dimension, 
which gives about a factor of 16 in computation time between 
the highest and lowest resolution. The results are very similar 
to the 2D case for the eccentricity and migration. The eccentric- 
ity damping time scale shortens slightly (by about 10 %) when 
comparing the highest and lowest resolution runs, while the mi- 
gration rate hardly changes at all. The inclination damping time 
shortens as well upon increasing the resolution (Fig. [T6l ), but 
this is largely a result of the change in eccentricity damping time 
(see below). Additional resolution studies performed in 3D for 
eo = 0.2, io - and eo = 0, io = 5 yielded very similar re- 
sults and are not shown here. In the case with eo - 0, io - 5 
degrees, we found that the inclination damping time was hardly 
affected by the changes in resolution, which is why the deviation 
in inclination damping times shown in fig. A. 2 is actually due to 
changes in the eccentricity damping time. 
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Fig. 16. Semi-major axis, eccentricity and inclination as a func- 
tion of time for a 20 niE planet in a three-dimensional disk, for 
different grid resolutions (see text). 



Both 2D and 3D resolution studies indicate that there are 
still residual changes at the 10 % level when the resolution 
changes, and this change mainly occurs in the eccentricity 
damping rate. Altogether our results demonstrate very clearly 
that the effect of an increased migration and inclination damp- 
ing for larger e and / in a locally isothermal disc is a ro- 
bust effect independent on resolution. The absolute magnitude 
will depend on physical details in the vicinity of the planet, 
where physics that we have not included in our models (e.g. 
thermal and radiative effects, MHD turbulence) will play a 
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role amongst others dNelson & Papaloizoull2004l: lNelsoDll2005l: 
iPaardekooper & Mellem"all2006l) 
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